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ABSTRACT 

An  efficient  method  is  presented  for  calculating  chemical  composition  in  a  perfectly 
stirred  reactor  at  a  given  pressure,  stoichiometry,  and  mass  flow  per  unit  volume.  Either 
temperature  or  enthalpy  may  be  prescribed  as  the  additional  condition  for  the 
computations.  By  employing  linearization  techniques,  the  nonlinear  equations  of  detailed 
finite  rate  chemical  kinetic  schemes  are  reduced  to  a  system  of  algebraic  equations  which 
are  solved  iteratively.  No  difficulties  are  experienced  in  obtaining  converged  solutions  using 
the  techniques  described  in  this  report.  A  computer  program  is  presented  for  solution 
of  problems  with  arbitrary  fuel/oxidant  combinations. 
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NOMENCLATURE 

Aj  Pre-exponential  coefficient  in  the  jth  reaction  rate  constant 

ak,2  The  kth  specie  coefficient  in  the  conservation  equation  for  the  Eth  specie 

Bj  Back  reaction  rate  of  the  jth  reaction 

bj  Back  reaction  rate  constant  in  the  jth  reaction 

Cg  Constant  in  the  Eth  conservation  equation 

Cp  Specific  heat  of  the  mixture 

CPi  Specific  heat  of  species  i 

dy  Third-body  efficiency  of  species  i  in  the  jth  reaction 

Ej  Activation  energy  in  the  jth  reaction  rate  constant 

Fj  Standard  molar  free  energy  of  species  i 

AFj  Change  in  free  energy  in  the  jth  reaction 

fj  Forward  reaction  rate  constant  in  the  jth  reaction 

H  Specific  enthalpy  of  the  mixture 

Hl  Enthalpy  lost  from  the  system 

hj  Specific  enthalpy  of  species  i 

Kj  Equilibrium  constant  for  the  jth  reaction 

KPj  Equilibrium  constant  for  the  jth  reaction  (pressure  units) 

m  Mass  flow  rate  into  the  reactor 

N  Number  of  linearly  independent  species 

NR  Number  of  reactions 

NT  Total  number  of  species 

P  Pressure 

R  Universal  gas  constant  (thermal  units) 
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R'  Universal  gas  constant  (pressure/volume  units) 

T  Temperature 

t  Time 

V  Volume  of  reactor 

W  Molecular  weight 

Xj  Third  body  in  a  dissociation  reaction 

Xjj  Coefficient  in  the  solution  matrix 

yj  Constant  in  the  solution  matrix 

a'|j,a'ij  Stoichiometric  coefficients  for  the  ith  species  in  the  jth  reaction 

j3j  Index  which  is  zero  for  second-order  and  unity  for  third-order  reactions 

71  Mole  mass  ratio  of  species  i  in  the  reactor 

71  Rate  of  production  of  species  i  due  to  chemical  reaction 

Sj  Finite  approximation  to  the  error  in  assumed  species  mole  mass  ratio 

6j  Finite  approximation  to  the  error  in  assumed  temperature  in  the  reactor 

p  Density  of  the  gas  mixture 

Tr  Residence  time  in  reactor  (average  value  for  all  molecules) 

SUPERSCRIPTS 

o  Denotes  standard  value  or  feed  quantity 

*  Denotes  approximate  value  of  a  variable 

SUBSCRIPTS 

i,k,C  Dummy  subscripts  on  species 

j  Denotes  reaction  number 

T  Denotes  temperature  coefficient 
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SECTION  I 
INTRODUCTION 

A  well-stirred  reactor  is  a  device  into  which  fuel  and  oxidant  are  injected,  mixed, 
and  reacted.  If  the  mixing  time  is  much  shorter  than  the  total  residence  time  in  the  reactor, 
then  the  temperatrue  and  composition  in  the  device  may  be  considered  perfectly  stirred, 
i.e.,  uniform  and  the  same  as  the  exhaust  from  the  device. 

The  concept  of  a  well-stirred  reactor  was  initially  developed  by  Longwell  and 
co-workers  (Refs.  1  through  3)  for  investigation  of  reaction  fuel/oxidant  systems.  Typically, 
measurements  were  made  of  injection  rates  that  caused  extinction  of  the  flame  for  various 
fuels,  equivalence  ratios,  pressures,  and  initial  temperatures.  These  measurements  were  used 
to  deduce  maximum  volumetric  rates  for  the  aggregate  of  the  participating  reactions,  i.e., 
the  "volumetric-heat-ielease  rates."  Their  lead  has  been  followed  by  several  groups  of 
workers,  namely:  Clarke,  et  al  (Refs.  4  through  6),  Hottel,  et  al  (Refs.  7  and  8).  A  large 
volume  of  literature  exists  in  the  chemical  process  industry  where  the  stirred  reactor  is 
used  extensively.  It  is  in  this  industry  where  stirred  reactor  theory  was  extended  to  more 
than  a  one-step  reaction  from  reactants  to  products  and,  in  particular,  incorporation  of 
statistical  models  of  turbulence  and  effects  of  incomplete  mixing/ stirring  (Refs.  9  through 
11).  Jenkins  and  Yumlu  (Ref.  12)  introduced  the  concept  of  making  detailed  reaction 
rate  studies  with  a  stirred  reactor.  Detailed  numerical  studies  were  carried  out  for  the 
hydrogen/ oxygen  system,  using  a  program  developed  by  Jones  (Ref.  13)  for  solving  the 
finite  rate  chemical  kinetic  equations.  The  computer  program  was  used  to  determine 
corrections  to  reaction  rate  constants  so  that  predicted  flame  temperatures  could  be  made 
to  agree  with  measured  values  for  a  limited  range  of  operating  conditions. 

The  strong  similarity  of  a  stirred  reactor  to  the  gas  turbine  engine  primary  combustor 
has  led  to  its  application  to  gaseous  pollutant  studies  (Refs.  14  through  16).  However, 
the  large  number  of  variables,  13  species  plus  temperature  even  for  methane  combustion, 
demanded  a  high  degree  of  sophistication  in  iterative  techniques  to  obtain  solutions.  Several 
solution  techniques  have  been  developed  (Refs.  13  and  16).  The  method  of  Ref.  13  is 
based  on  solving  a  set  of  nonlinear  algebraic  equations  which  describe  the  chemical  kinetics 
of  a  fuel/oxidant  system.  The  equations  are  linearized  about  the  solution  and  an  iterative 
procedure  used  to  obtain  converged  solutions.  A  similar  technique  is  described  in  this 
report  including  extension  to  variable  temperature  as  well  as  specie  composition.  The 
method  used  in  Ref.  16  is  referred  to  as  a  technique  for  solving  the  time  dependent 
form  of  the  stirred  reactor  equations  by  a  relaxation  method.  Another  method  employs 
numerical  integration  techniques  with  the  time  dependent  term  in  the  stirred  reactor 
equations  relaxed  to  zero  when  the  solution  is  obtained.  Considerable  difficulty  has  been 
experienced  with  the  iterative  method  of  solving  the  linearized  algebraic  equations.  The 
equations  are  highly  nonlinear  and  the  generally  used  Newton-Raphson  technique  requires 
an  approximate  solution  from  which  to  start  the  iteration.  The  difficulty  increases  rapidly 
as  the  number  of  variables  increases.  The  convergence  problem  becomes  quite  severe  in 
chemical  laser  calculations  where  each  excited  state  of  a  specie  is  considered  as  a  separate 
specie  (Ref.  17).  The  perfectly  stirred  reactor  may  constitute  an  ideal  chemical  laser  because 
of  the  uniformity  of  chemical  state  over  the  whole  volume  and  its  ability  to  operate 
steady  state  as  far  from  chemical  equilibrium  as  desired  (providing  flameout  does  not 
occur). 
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The  technique  and  program  presented  in  this  report  have  been  used  for  a  variety 
of  problems  over  a  large  range  of  input  pressures,  temperatures,  stoichiometries,  and  mass 
flows/unit  volume.  Convergence  is  obtained  without  difficulty  even  for  large  numbers  of 
variables,  and  parametric  calculations  have  been  carried  out  with  very  short  computational 
time  requirements. 

The  system  of  equations  to  be  solved  and  the  method  of  solution  is  described.  A 
listing  of  the  program  and  a  description  is  then  given  with  representative  sample  calculations 
for  a  hydrogen/air  and  a  methane/air  system. 

SECTION  II 

THEORY  OF  THE  IDEAL  STIRRED  REACTOR 
WITH  FINITE  RATE  CHEMISTRY 

Suppose  that  a  reactant  gas  stream  flows  at  a  mass  flow  rate  (m)  into  a  reactor 
of  volume  V,  its  composition  being  defined  by  a  set  of  values  (7I),  representing  the  inlet 
concentration  of  the  component  species  i.  The  composition  in  the  reactor  is  assumed 
uniform  (perfectly  stirred)  and  identical  with  that  in  the  exhaust  stream;  it  is  defined 
by  the  set  of  values  (7i).  The  temperature  T  is  also  supposed  uniform  in  the  reactor 
and  exhaust  stream.  The  reactor  walls  are  assumed  to  be  entirely  inert  chemically  and 
do  not  catalyze  the  reactions  in  any  way.  The  gases  are  assumed  to  be  thermally  perfect. 

Conservation  equations  can  be  written  as  follows: 

Species  i 


-y{)  +  Yi  - 


d>'i 

dT 


The  species  production  rate  (71)  is  obtained  from 


NR 

v.  =  X  (a.'.'-a'XF.-  B;) 
'  1  j_]  ij  >r  )  y 


(1) 


(2) 


Energy 


NT 


dH 


^£W-yW-ir  +  HL  O) 

For  steady-state,  adiabatic  operation,  these  equations  reduce  to  those  of  Ref.  12: 


p-(yi0-7i)  +  *1-  0 


(4) 


hl  =  0 
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The  ay'  and  ay  are  stoichiometric  coefficients  in  the  jth  chemical  reaction,  Fj  is  the 
rate  of  the  forward  reaction,  and  Bj  is  the  rate  of  the  backward  reaction.  The  quantities 
hj  and  hi  are  specific  enthalpies  (functions  of  temperature  only)  in  the  reactor  and  inlet 
stream,  respectively,  of  species  i.  The  steady-state  form  of  the  equations  is  used  in  this 
report. 


The  reaction  rates  are  related  to  the  concentrations  in  the  reactor  by 


Fi  " 


fixr 


aij 


and  B. 


/8,  NT 

bixi  ‘I,w 


(5) 


where  p  is  the  density  of  the  gas  mixture  in  the  reactor,  fj  and  bj  are  rate  constants 
of  the  forward  and  backward  reactions,  respectively,  and  Xj  is  a  third  body  in  a  dissociation 
reaction  and  is  given  by 


Xi 


NT 

£  d-.y 

'Y' 


(6) 


where  the  dy  represent  third-body  efficiencies.  The  rate  constants  are  of  the  form 

fj  -  A-t"*  exp  (-E./RT)  (7) 

however,  the  back  reaction  rate  constant  is  not  chosen  independently.  The  rate  constants 
obey  the  following  relation: 


fjKj 

(8) 

NT 

Kp  (R"R  1 

(9) 

exp(-AFj/RT) 

(10) 

r  nt  i 

i°/RTl 

(11) 

Kj  is  the  equilibrium  constant  and  F?  is  the  molar  free  energy  of  species  i.  The  system 
of  equations  is  completed  with  the  equation  of  state 


NT 

P-PjT-  pR'T.2  y. 


(12) 


where  P  is  the  system  pressure  and  W  is  the  molecular  weight  of  the  gas  mixture. 

The  equations  form  a  complete  algebraic  system  of  nonlinear  equations  for  steady-state 
operating  conditions  and  can  be  solved  by  standard  techniques  once  a  reaction  kinetic 
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scheme,  rate  data,  and  thermodynamic  data  are  given.  An  efficient  method  of  solution 
is  described  which  has  proved  quite  successful  for  a  large  range  of  problems. 


It  is  assumed  that  at  least  one  approximate  solution  can  be  obtained  using  an 
equilibrium  composition  calculation,  a  one-dimensional  finite  rate  chemistry  calculation 
(plug  flow)  or  solution  by  integration  of  the  transient  stirred  reactor  Eq.  (1)  with 
the  transient  term  relaxed  to  zero.  Guessing  a  solution  is  nearly  impossible  for  all  but 
very  simple  reaction  schemes  such  as  that  of  Ref.  12. 


These  approximate  values  (7?)  and  T*  are  used  in  a  Newton-Raphson  procedure  to 
obtain  first-order  correction  terms  from  which  new  values  of  7*  and  T  can  be  obtained 


y *  +  si 

(13) 

* 

T  + 

(14) 

By  repeating  the  procedure,  the  6;  and  5f  can  be  relaxed  to  some  negligible  value;  thus 
final  values  (solutions)  for  7*  and  T  are  obtained. 


The  correction  equations  are  obtained  as  follows: 


y.  =  y.*  +  dy.  and  Yi  “  y  *  + 


N  T  dy- 


dy.  -  irr  dT  +  dP  +  2  ^  dx  k 

<5t  dp  ^  k=i  dyk 


?!i 

i  dT 


dF; 


NT  dF. 


dF:  =  -*JL  dT  +  — L  ■  v’ 


, -  dp  +  2  3 - 

op  k=l  dyk 


dB, 


dB. 


NTdB. 


-  FT dT  ♦  or  * +  Am 

NT  dX:* 

dX.  =  X  3—  dyk 
J  k=l  dyk 


dyk 


dyk 


Not  all  variables  are  linearly  independent;  thus, 

dp  =  dT  +  V 

y  dT  k=i  dyk  rk 


dyg  =  2  ak  gdy^  £  =  N  +  1,  NT 


y?  =  C£  +  ak,£Fk 


(15) 

(16) 

(17) 

(18) 

(19) 

(20) 

(21) 

(22) 


The  Cg  values  are  obtained  from  element  conservation  equations. 
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The  partial  derivatives  are  obtained  from  the  preceding  conservation  equations.  If 
temperature  is  an  unknown  variable,  the  energy  equation  is  used  to  obtain  a  correction 
equation  for  it: 


H 

NT 
-  £, 

(23) 

dH 

NT  N  /  NT  \ 

■  *  £i(k‘  ♦  £,+1  •k*)^ 

(24) 

N  /  NT  \ 

‘  CPJT  +  +  £„♦!  '*■*• « )*' 

(25) 

where 

CP 

NT 

-  £  y-Cn 
i=l  r'  Pi 

(26) 

Summarizing  gives 

y.  =  y;(yk,  D  i,k  =  1,N 

(27) 

P  =  p(yk ,  T) 

(28) 

yg  -  yg(yk,T)  t  =  n+i,nt 

(29) 

yi°-yi*  +  rR*yi*  -  -  y^ 

(30) 

where 

pv  •  j 

rR  *  —  “  rR  +  drR 
m 

(31) 

and 

jp 

rR  ~  P 

(32) 

Substituting  for  drR  gives 


y  t  ~  y i  +  rR  Yi 


-  ^  -  'H^i  ♦  +  Y 


(33) 


Substituting  for  d-ft  and  d7j 


O  *  •  * 

yj  -  Yi  +  rR  yi 


"rR* 


*y5 

(?T 


-r 


R 


-  yV(l  * 


NT 

E=N+1 


(34) 
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where 


and 


2.1  METHOD  OF  SOLUTION 


(35) 


When  a  converged  solution  is  obtained  both  the  left-  and  right-hand  sides  of  Eq. 
(34)  are  approximately  equal  to  zero.  By  using  the  Newton-Raphson  technique,  estimates 
(7i)  for  N  of  the  total  number  of  species  NT  are  required  together  with  an  estimate  of 
the  temperature  (T*).  These  estimates  are  used  to  evaluate  the  terms  on  the  right-hand 
side  of  Eq.  (34).  Thus,  a  matrix  can  be  formed  as  follows: 


Yi 


N 

-  £ 
k=i 


xi,ksk 


+  x.  i  —  1,N 


(36) 


where  dy*  is  replaced  by  the  finite  approximation  (St),  dT  is  replaced  by  St,  yi  represents 
the  left-hand  side  of  Eq.  (34),  and  the  xik  represent  the  coefficients  on  the  right-hand 
side.  The  (N  +  1)  row  of  coefficients  to  be  added  to  the  above  matrix  are: 


yN+1  =dH=AH=H0-H 

(37) 

NT 

N+l,k  =  hk  +  ^I+1  flk,£h£ 

(38) 

yN+1  =  k=l  XN+l.k^k  +  CpST 

If  the  stirred  reactor  composition  is  required  at  a  specified  temperature  and  pressure, 
instead  of  in  terms  of  the  inlet  energy  and  pressure,  then  only  the  coefficients  of  an 
N  x  N  matrix  are  required,  and  the  St  terms  are  ignored. 

The  matrix  is  solved  for  the  unknown  6*  and  St,  and  values  of  Tk  -  7k  +  and 
T  =  T*  +  ST  are  computed.  The  process  is  repeated  until  6k  «  7k  and  St  «  T. 

Some  difficulty  was  experienced  in  obtaining  converged  solutions  for  N  >  6.  The 
primary  problem  was  found  to  be  the  exponential  response  of  specie  derivative  terms 
when  large  temperature  corrections  were  required  in  constant  enthalpy  calculations.  This 
problem  was  overcome  by  deleting  all  temperature  coefficients  from  the  matrix  and 
allowing  several  composition  iterations  to  be  made  between  each  temperature  correction 
obtained  as  follows: 

Let  H0  -  H  =  CpAT 
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H  and  Cp  are  evaluated  from  the  current  value  of  y\  and  T*,  and  a  new  value  of  T 
is  given  by 


T  =  T*  +  AT  =  T*  +  (H0  -  H)/C 

It  was  also  found  that  if  an  odd  fraction  was  used  to  multiply  each  indicated  correction 
for  species  and  a  large  fraction  used  to  reduce  the  temperature  correction  until  the  AT 
changed  sign  and  then  use  a  small  fraction  to  reduce  it,  oscillations  about  the  solution 
were  eliminated. 

2.2  REACTION  RATE  SENSITIVITY  COEFFICIENTS 

The  experiments  described  in  Ref.  12  were  designed  to  furnish  data  to  allow 
adjustment  of  the  reaction  rates  so  that  the  reactor  temperature  predicted  by  the  computer 
program  could  be  made  to  agree  with  the  measured  value.  An  analysis  of  the  sensitivity 
of  predicted  reactor  temperature  to  individual  reaction  rate  constants  was  carried  out  to 
determine  which,  if  any,  had  a  dominant  role  in  determining  the  reactor  temperature. 
Since  reactions  rates  are  functions  of  composition  also,  the  calculations  were  repeated 
with  different  values  of  initial  water  concentration  in  the  feed  of  hydrogen/oxygen.  The 
sensitivity  coefficient  is  defined  as 


with  units  of  °K.  It  represents  the  change  in  reactor  temperature  for  an  infinitesimal  change 
in  reaction  rate  constant  (fj).  The  quantity  is  evaluated  numerically  for  an  0.1-percent 
increase  in  reaction  rate  constant  (fj): 


dr 

i 

i 


f.*I 

Af. 

] 


This  calculation  is  carried  out  for  each  reaction  rate  in  sequence.  The  sensitivity  to 
initial  feed  composition  may  also  be  evaluated  with  the  computer  program.  Control  integers 
are  incorporated  whereby  the  feed  composition  of  any  variable  is  assigned  multiples  of 
5  percent  of  the  inert  species  composition.  Reaction  rate  sensitivity  coefficients  are 
re-evaluated  for  each  new  feed  composition.  The  range  of  possible  experiments  is 
considerably  increased  by  this  means  since,  for  the  methane/air  system,  CO,  CO2,  and 
H2O  can  be  readily  included  in  the  feed  stream. 


This  subroutine  is  used  to  gain  insight  into  which  reactions  are  rate  controlling  for 
various  reactor  temperatures  and  feed  compositions.  Tabulating  the  sensitivity  coefficients 
for  a  wide  range  of  input  conditions  and  reactor  residence  times  shows  whether  some 
reactions  may  even  be  eliminated  in  practical  calculations.  If  the  sensitivity  coefficient 
in  a  reaction  is  always  small  relative  to  others,  then  this  reaction  could  probably  be 
eliminated.  An  example  is  discussed  in  a  later  section  for  the  methane/air  system.  An 
example  for  the  hydrogen/ oxygen/steam  system  is  given  in  Ref.  12. 
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SECTION  III 
COMPUTER  PROGRAM 

A  computer  program  was  written  for  the  Raytheon  520  digital  computer  in  Fortran 
IV  language.  A  listing  of  the  program  is  given  in  Appendix  III.  The  program  is  written 
in  -a  general  format  so  that  any  chemical  system  can  be  investigated.  The  maximum  number 
of  different  species  is  20,  and  the  maximum  number  of  chemical  reactions  is  76.  Either 
of  these  numbers  can  be  readily  increased  or  decreased  as  desired. 

A  number  of  options  are  included  which  allow  a  range  of  known  flame  temperatures 
and/or  a  range  of  mass  flows  per  unit  volume  to  be  used.  An  additional  subroutine  has 
been  incorporated  to  calculate  temperature  sensitivity  coefficients  for  reaction  kinetic 
schemes. 

The  program  is  divided  into  several  subroutines  which  execute  specific  parts  of  the 
overall  calculation.  The  calculations  are  controlled  by  integer  variables  defined  in  the  main 
program.  Some  comments  have  been  added  to  help  identify  the  types  of  calculations  carried 
out  in  the  program. 

3.1  MAIN  PROGRAM 

Data  cards  are  as  follows: 

I.  P,  T.  TO,  DELT,  TM1N 

II.  MDOTV,  DMDV,  NT,  N,  LT, 

MT,  IT,  NR,  NN,  NNT 

III.  XO(l),  I  =  1,  NT 

IV.  X(l),  I  =  1,  NT 

V.  W(l),  I  =  1,  NT 

I.  P  is  the  system  pressure  (atm),  T  is  the  known  or  estimated  flame  temperature 
(°K),  TO  is  the  feed  temperature,  DELT  is  the  temperature  decrement  desired  when 
temperature  is  not  a  variable,  and  TMIN  is  the  minimum  temperature  at  which  a  calculation 
is  to  be  attempted. 

II.  MDOTV  is  the  mass  flow  per  unit  volume  (g  cm’3  sec’1 ),  DMDV  is  the  factor 
by  which  MDOTV  is  to  be  increased  (or  decreased)  in  parametric  calculations,  NT  is  the 
total  number  of  chemical  species,  N  is  the  number  of  linearly  independent  species,  LT 
is  the  number  of  times  MDOTV  is  to  be  multiplied  by  DMDV,  MT  is  the  number  of 
times  DELT  is  to  be  subtracted  from  T,  IT  is  a  control  variable  (=1  for  constant 
temperature,  and  =2  for  constant  enthalpy  calculations).  NR  is  the  total  number  of  chemical 
reactions,  NN  is  a  control  variable  (=0  to  bypass  the  reaction  rate  sensitivity  calculation 
and  =1,  the  species  number  whose  concentration  is  to  be  changed  in  the  reaction  rate 
sensitivity  calculations),  and  NNT  is  the  number  of  times  species  NN  is  to  be  changed. 
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III.  XO(I)  is  the  molar  ratio  of  species  I  in  the  feed.  The  order  of  the  species 
is  not  arbitrary;  the  species  to  be  evaluated  with  conservation  equations  must  be  last. 

IV.  X(I)  is  the  molar  ratio  of  species  I  in  the  reactor.  The  order  of  the  numbers 
must  be  the  same  as  in  the  card  above.  This  requirement  applies  to  all  specie  and  specie 
data  cards. 

V.  W(I)  is  the  gram-molecular  weight  of  species  I. 

Calculations  carried  out  in  the  main  program  are: 

i.  normalization  of  input  concentrations  followed  by  conversion  to  mole 
mass  units  and  summed  tc  obtain  the  reciprocal  molecular  weight  of 
the  mixture. 

ii.  Subroutine  BB3  is  called.  In  this  subroutine,  species  enthalpies,  specific 
heats,  free  energies,  and  derivatives  of  free  energies  are  evaluated.  The 
first  time  BB3  is  called,  the  initial  enthalpy  corresponding  to  the  feed 
temperature  and  composition  is  calculated  (HO).  The  next  time  BB3 
is  called,  the  current  value  of  T  is  used.  The  subroutine  is  not  called 
again  until  T  is  changed.  The  free  energies  and  their  temperature 
derivatives  are  used  to  evaluate  equilibrium  constants  and  their 
temperature  derivatives. 

iii.  Calculation  of  feed  and  estimated  flame  densities. 

iv.  Calculation  of  estimated  reactor  residence  time. 

3.2  SUBROUTINE  COMP 

This  subroutine  calls  a  number  of  subroutines  in  which  the  data  necessary  to  evaluate 
the  coefficients  for  the  iteration  matrix  are  calculated.  When  the  coefficients  have  been 
evaluated,  a  matrix  reduction  routine  is  called  which  returns  the  values  of  the  corrections 
to  the  initial  estimates.  A  check  is  then  made  to  determine  whether  the  largest  correction 
would  produce  a  negative  species  concentration,  if  so  then  all  correction  terms  are  reduced 
by  the  amount  necessary  to  prevent  this,  and  finally  the  corrections  are  multiplied  by 
0.37.  This  has  ensured  satisfactory  convergence  for  all  calculations  attempted.  The  new 
values  of  the  species  are  evaluated,  and  a  check  is  made  to  see  if  the  original  equations 
are  satisfied;  if  so,  the  converged  solutions  are  used  to  compute  final  values  of  molecular 
weight,  density,  and  enthalpy  and  these  values  are  returned  to  the  main  program. 

The  sequence  of  calculation  is  as  follows: 

i.  Subroutine  BB1  is  called  once  only,  to  obtain  coefficients  for  the 
conservation  equations  and  their  derivatives. 

ii.  Subroutine  BB4  is  called  every  time  a  new  set  of  composition 
corrections  is  to  be  calculated.  This  subroutine  also  calls  Subroutine 
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RATES  if  the  temperature  is  changed,  otherwise  once  only  to  establish 
the  values  of  the  reaction  rate  constants  and  their  derivatives.  The 
majority  of  the  partial  derivatives  are  calculated  in  subroutine  BB4. 

iii.  The  coefficients  for  the  matrix  are  evaluated. 

iv.  Subroutine  CHOLES,  the  matrix  reduction  routine,  is  called  to  evaluate 
the  values  of  the  corrections. 

v.  The  corrections  are  used  and  convergence  to  the  solution  evaluated. 

a3  SUBROUTINE  BB1 

This  subroutine  evaluates  the  coefficients  of  the  conservation  equations  using  constants 
which  are  read  in  the  first  time  the  subroutine  is  called.  The  data  are  read  in  as  follows: 

I.  J(I),  I  =  NL  where  NL  is  the  number  of  linearly  dependent  species  (NL  =  NT 
-  N).  Each  integer  read  in  denotes  how  many  coefficients  are  to  be  called  for  in  the 
J(I)  th  conservation  equation. 

II.  (NS(K,I),  SC(K,  I),  K  =  1,  JJ),  I  =  NL  where  JJ  (=J(I))  is  the  number  of  species 
in  each  conservation  equation.  NS  (K,  I)  is  the  number  of  the  species,  and  SC(K,  I)  is 
the  value  of  the  coefficient  in  the  conservation  equation. 

An  example  for  conservation  of  H2  and  O2  in  the  H2/O2  system  is 


h2  = 

1/2  H  + 

1/2  OH  +  H20  +  H2 

02  = 

1/2  0  + 

1/2  OH  +  1/2  H20  -t 

■  02 

Card  I 

4,  4 

Card  Ha 

1,  0.5 

O 

wo 

O 

•V 

CO 

5,  1.0 

Card  lib 

2,  0.5 

3,  0.5  4,  0.5 

6,  1.0 

where  1  =  H,  2  =  O,  3  =  OH,  4  =  H2O,  5  =  H2,  6  =  O2.  The  order  of  the  cards  must 
correspond  to  the  order  of  the  molecular  weights  in  the  main  program. 

3.4  SUBROUTINE  BB3 

This  subroutine  is  used  to  calculate  species  enthalpy,  specific  heat,  free-energy,  and 
derivatives  of  free  energy.  The  first  time  subroutine  BB3  is  called,  a  set  of  coefficients 
is  read  in  as  follows: 

C(I,  KK),  I  =  1,  14  KK  =  1,  NT 
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These  coefficients  correspond  to  those  given  in  NASA-SP-3001  (Ref.  18).  The  first  seven 
coefficients  correspond  to  the  temperature  range  from  300  to  1000°K  and  the  second 
seven  are  for  the  range  from  1000  to  5000°K,  i.e.,  fourteen  coefficients  per  specie.  These 
data  cards  must  correspond  to  the  order  of  the  molecular  weights  in  the  main  program. 

3.5  SUBROUTINE  CHOLES 

This  is  a  standard  matrix  reduction  routine  developed  by  Glassman1.  It  is  based  on 
the  Choleski  method  for  matrix  factorization. 


&6  SUBROUTINE  BB4 


This  subroutine  contains  the  "chemical  kinetics"  of  the  chemical  system  being 
investigated.  The  data  cards  contain  all  the  information  relating  to  chemical  reactions  and 
reaction  rate  constants.  Subroutine  RATES  is  called  first  to  evaluate  reaction  rates  and 
reaction  rate  derivatives;  then  the  partial  derivatives  required  for  the  matrix  reduction 
equations  are  evaluated. 

The  following  nomenclature  is  used  to  denote  the  corresponding  variables  in  Section 
III. 

dVi 

BO(I)  =  y.°  B(I)  s  yj  R  =  p  DBDBU.K)  = 

dy •  dy • 

DBDT(I)  ^  ^  DBDR(I)  ^  BD(D  3  y. 

The  reaction  rate  cards  are  read  in  as  follows: 

I.  NR 

II.  ITR(IR),  ITP(IR),  (KK(IR,  I),  I  =  1,  4),  A(IR), 

EX(IR),  TEX(IR),  IE(IR),  IC(IR),  IB(IR) 

III.  (I KIR),  ER(IR,  I),  I  =  1,  4) 

IV.  KS(IR,  I),  I  =  1,  IC  (IR) 


I.  NR  is  the  total  number  of  reactions,  IR  is  an  integer  in  the  range  IR  =  1,  NR. 

II.  ITR(IR)  and  ITP(IR)  are  the  number  of  reactants  and  products  in  the  IRth 
reaction.  KK(IR,  I),  I  =  1,  4  denotes  the  Ith  species  in  the  IRth  reaction,  and  a  maximum 
of  four  is  allowed  (in  a  recombination  reaction,  e.g.,  H  +  H  +  M  =  H2+M,  M  has 
the  value  NT  +  1  and  the  four  integers  corresponding  to  the  reactants  and  products  1, 
1,  5,  7  in  the  nomenclature  of  Subroutine  BB1).  A(IR)  is  the  pre-exponential  factor, 
EX(IR)  is  an  activation  energy,  and  TEX(IR)  is  a  temperature  exponent  in  the  IRth  reaction 
rate  constant. 

1H.  N.  Glassman,  Central  Computer  Operations,  ARO,  Inc.,  AEDC 
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IE(IR),  IC(IR),  and  IB(1R)  are  integer  controls  for  different  types  of  reactions.  They 
are  most  easily  explained  by  the  following  examples: 

a.  Two  body  reaction  with  forward  rate  constant  data  read  in:2 

H  +  02  =  OH  +  O 

1R  =  3,  A(l)  =  1  9  x  10“  EX(1)  =  9007.0,  TEX(l)  =  0 
IE(1)  =  0,  IC(1)  =  0,  IB(1)  =  1 
The  data  card  for  this  reaction  is 

2,  2,  1,  6,  3,  2,  1.9E  +  14,  9007.,  0.,  0,  0,  1 
Cards  III  and  IV  are  not  required. 

b.  Three-body  reactions  must  be  written  as  recombination  reactions: 

H+H  +  M  =  H2+M 

If  the  known  rate  is  kb  for  the  back  reaction,  then  IB(1)  =  -1  so  that 
the  reaction  rate  corresponds  to  the  right-hand  side. 

The  card  would  then  read 

3,  2,  1,  1,  4,  7,  5.  0E  +  18,  0.0,  -1.15,  0,  0,  -1 

c.  Three-body  reaction  with  H20  4.7  times  as  efficient  as  any  other  species 

acting  as  a  third  body:  IR  =  8  H  +  OH  +  M  =  HzO  +  M 

Card  II  would  read  (with  IE(2)  =  1) 

3,  2,  1,  3,  4,  7,  1.2E  +  17,  -252.,  0.,  1,  0,  1 
Card  III  would  read 

4,  3.7  where  11(2)  =  4  and  ER(2)  =  3.7 

The  value  of  ER(IR)  is  always  1.0  less  than  the  third-body  efficiency.  If 
they  are  all  equal  to  unity,  the  card  is  not  required.  If  the  reverse  reaction 
is  effectively  zero,  then  IB(2)  =  0,  and 

Card  II  would  read 

3,  2,  1,  3,  4,  7,  1.2E  +  17,  -252.,  0.,  1,  0,  0 

d.  Special  Reaction: 

IR  =  9  H  +  02  +  M  -*■  H  02  +  M  2  OH  +  M 
This  reaction  is  handled  with  data  given  in  Card  IV. 

Reactions  and  late  data  are  detailed  in  Table  1  (Appendix  II). 
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II.  3,  2,  1,  6,  7,  2. IE  +  18,  0.,  -1,  1,  5,  0 

III.  4,  2.6 

IV.  1,  6,  7,  3,  3 

where  IC(3)  =  5,  the  number  of  molecules  involved  whose  numbers  are 
specified  in  Card  IV.  Card  III  is  required  since  H2O  is  assumed  to  be  3.6 
times  as  efficient  as  all  other  molecules  acting  as  third  bodies. 

3.7  SUBROUTINE  RATES 

This  subroutine  evaluates  the  values  of  the  reaction  rate  constants  and  derivatives 
required  in  subroutine  BB4. 

In  Section  III,  fj  =  Aj  Tnj  exp  (-  Ej/RT) 

in  this  subroutine  fj  =  FR(IR)  =  A(IR)  Ttex(ir7  eE  <IR  >/T 

In  Section  III,  bj  =  fjKj  (Kj  =  EQ  in  this  subroutine). 

This  equilibrium  constant  EQ  is  evaluated  using  free  energies  obtained  in  Subroutine  BB3. 

3.8  SUBROUTINE  RCHEK 

Reaction  rate  sensitivity  coefficients  are  obtained  in  this  subroutine  using  a  numerical 
approximation  for 


The  subroutine  is  called  if  the  integer  NN  f-  0  in  the  main  program.  If  NN  is  nonzero, 
it  is  equal  to  the  number  of  one  of  the  constituent  species. 

The  control  integer  NNT  determines  if  and  how  many  times  the  feed  concentration 
of  species  NN  is  to  be  increased  as  a  fixed  percentage  (5  percent  of  the  last  species 
concentration  with  the  inert  N2  the  species  NT).  This  enables  the  sensitivity  coefficients 
to  be  determined  as  a  function  of  each  species  thought  to  be  significant.  Inclusion  of 
the  call  statement  for  this  subroutine  within  the  loops  for  temperature  and/or  mass  flows 
per  unit  volume  enables  the  effect  of  temperature  and  residence  time  to  be  calculated 
also. 


The  sensitivity  coefficient  is  evaluated  as  follows: 

An  integer  (ICR)  is  specified  and  when  subroutine  rates  is  called,  the  pre-exponential 
coefficient  A(IR)  is  relabelled  A(ICR)  and  multiplied  by  1.001.  The  calculation  is  carried 
out  with  IT  =  2,  i.e.,  at  constant  enthalpy  and  a  new  value  of  temperature  determined 
and  compared  with  the  value  calculated  previously  (before  A(IR)  was  multiplied  by  1.001). 
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The  difference  in  these  two  temperatures  is  multiplied  by  1000  and  this  represents  the 
temperature  change  AT  for  a  0. 1-percent  change  in  reaction  rate  constant.  The  sequence 
is  carried  out  for  ICR  =  1,  (NR  +  1),  i.e.,  all  the  reactions  are  investigated  including 
the  (NR  +  l)th,  which  is  a  check  to  see  if  AT  =  0.  for  this  case. 

This  subroutine  is  useful  for  determining  if  certain  reactions  have  a  negligible  effect 
on  the  chemical  system  (AT  «  0  for  a  range  of  stoichiometries,  residence  times, 
temperatures,  etc.).  If  they  are  negligible,  they  can  be  eliminated  to  speed  up  the  calculation 
and  reduce  computational  time. 

3.9  SAMPLE  OUTPUT 

The  sample  output  data  immediately  after  the  listing  are  for  a  CH4  /Air  system  detailed 
in  the  next  section.  The  integer  controls  correspond  to  the  use  of  subroutine  RCHEK. 
Thermodynamic  data  evaluated  at  temperature  T0  °K  are  listed  on  the  first  page  of  output 
lists.  The  coefficients  used  in  Subroutines  BB1  and  BB4  are  listed  on  the  second  page. 
Solutions  for  the  case  7CH3  (BO(6))  equals  0.0  and  0.5  percent  of  7ch4  (BO(5))  are 
presented  on  the  third  page.  The  first  solution  is  an  IT  =  1  case  to  determine  the  enthalpy 
at  the  given  temperature  and  is  written  from  the  main  program.  The  second  solution  (IT 
=  2)  contains  the  sensitivity  coefficients  with  the  same  feed  as  for  the  first  solution  and 
is  written  from  Subroutine  RCHEK.  The  third  solution  (IT  =  2)  is  written  from 
Subroutine  RCHEK  for  7c h 3  =  0.005  7ch4. 

SECTION  IV 

SAMPLE  CALCULATIONS  AND  DISCUSSION 
4.1  HYDROGEN/OXYGEN  SYSTEM 

The  chemical  kinetic  scheme  and  reaction  rate  data  used  are  given  in  Table  I.  Enthalpy 
versus  temperature  data  are  shown  in  Figs,  la  and  b  (Appendix  1)  for  equivalence  ratios 
(ER)  =  0.1  and  1.0.  The  minimum  in  each  curve  corresponds  to  the  temperature  at  which 
the  flame  "blows  out."  The  stirred  reactor  requires  a  significant  rate  of  production  of 
the  "active"  species  H,  O,  OH,  and  H2O  to  sustain  a  flame,  and  a  minimum  residence 
time  is  required  for  this  to  occur  in  the  reactor.  This  is,  of  course,  why  the  minimum 
shifts  to  higher  temperatures  as  the  residence  time  is  decreased  (increased  m/v). 

The  chemistry  in  perfectly  stirred  reactors  (PSR)  differs  somewhat  from  that  in 
premixed  plug  flow  reactors  (PFR).  In  a  PFR,  the  H,  O,  and  OH  species  concentrations 
increase  exponentially  with  time  during  the  ignition  process,  followed  by  a  decrease  as 
the  temperature  increases  to  the  flame  temperature.  In  a  PSR,  the  production  rate  of 
H,  O,  and  OH  is  always  positive  as  required  by  the  species  conservation  equation.  This 
requirement  leads  to  different  behavior  with  increasing  time  as  shown  in  Figs.  2a  through 
c.  This  characteristic  may  be  used  to  advantage  in  a  chemical  laser  where  a  supply  of 
radicals  is  required  in  a  chain-type  reaction,  operating  marginally  past  the  spontaneous 
ignition  point  (Ref.  17).  The  concentrations  as  a  function  of  temperature  with  enthalpy 
constant  in  both  reactor  models  are  essentially  identical  as  shown  in  Fig.  3  with  numerical 
errors  accounting  for  the  difference.  The  PFR  calculations  were  obtained  with  a  separate 
computer  program. 
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4.2  METHANE/AIR  SYSTEM 

The  chemical  kinetic  and  reaction  rate  data  are  given  in  Table  II.  Nitric  oxide 
formation  represents  a  composite  of  the  following  reactions: 

a.  N  +  02  =  NO  +  0  ba  =  1.55  x  109  T  e"19-450/7 

b.  0  +  N2  =  NO  +  N  fb  =  1.36  x  io 14  e"37.750/T 

c.  N2  +  02  =  NO  +  NO  fc  =  9.1  x  io24  T-2-5  e“64-300'/T 

The  rate  of  production  of  NO  is  given  by 

^2  *  =  2  *  =  *b>  ••c‘’  ^NO  =  2fbyN/02 

The  error  in  NO  concentration  is  small  for  temperatures  <  2500.0°K  and  could  be  readily 
reduced  with  some  numerical  experimentation  on  the  reverse  reaction  rate  constant. 

An  enthalpy  temperature  map  for  this  system  is  shown  in  Fig.  4.  Some  results  obtained 
using  the  reaction  rate  check  subroutine  are  shown  in  Table  III. 

The  results  are  presented  for  a  few  calculations  only.  In  practice,  a  comprehensive 
series  of  tables  would  be  necessary  with  as  large  a  range  as  possible  of  input  conditions, 
composition,  stoichiometry,  and  residence  time  (tr).  Inspection  of  the  table  shows  that 
the  nitric  oxide  production  equation  could  probably  be  eliminated  from  thermal  balance 
calculations  for  this  range  of  conditions  since  its  sensitivity  coefficient  is  always  relatively 
small.  It  obviously  should  be  retained  if  nitric  oxide  composition  is  to  be  estimated.  A 
marked  sensitivity  of  nearly  all  reactions  is  evident  near  the  blowout  condition  (see  Table 
III)  for  m/V  =  16,  T  =  1500°K.  This  table  does  not  indicate  any  operating  condition 
for  which  only  one  reaction  has  a  high  sensitivity  coefficient.  This  is  what  one  should 
look  for  in  a  comprehensive  series  of  tables  since  this  indicates  that  the  reaction  rate 
constant  of  this  reaction  can  be  adjusted  to  make  predicted  and  experimental  temperatures 
agree.  Repeating  the  procedure  for  each  experiment  with  one  or  more  dominant  reactions 
allows  successive  improvement  of  each  rate  constant  to  be  made. 

It  should  be  remembered  that  the  stirred  reactor  concept  is  limited  to  evaluating 
a  given  chemical  kinetic  and  reaction  rate  scheme,  i.e.,  it  cannot  serve  as  a  source  of 
fundamental  data  on  individual  chemical  reactions.  It  can,  however,  be  used  (numerically) 
to  evaluate  the  possible  influence  of  proposed  reactions.  Thus,  the  stirred  reactor  provides 
a  simple  means  of  evaluating  simplified  schemes  designed  to  approximate  more  detailed 
chemical  kinetic  schemes. 


SECTION  V 
CONCLUSIONS 

The  computer  program  has  been  used  successfully  to  solve  large  sets  of  nonlinear 
equations  in  an  efficient  manner.  No  difficulties  in  obtaining  solutions  are  experienced 
if  the  following  directions  are  followed  closely: 
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1.  Obtain  one  set  of  concentrations  from  a  finite  rate  solution  or  an 
equilibrium  solution. 

2.  If  temperature  is  to  be  varied,  restrict  the  temperature  change  to  less  than 
300°K.  If  larger  temperature  changes  are  desired,  then  the  solution  should 
be  "marched"  i.e.,  split  up  the  desired  change  into  a  number  of  smaller 
changes. 

The  program  allows  the  user  to  evaluate  a  chemical  system  parametrically  from 
flameout  to  equilibrium  with  any  desired  pressure,  temperature  or  enthalpy,  feed 
composition,  and  reactor  residence  time  via  the  mass  flow  rate  per  unit  volume.  The 
sensitivity  of  the  chemical  system  to  individual  reaction  rate  constants  and  feed  composition 
can  also  be  evaluated  for  prescribed  pressure,  temperature  or  enthalpy,  and  residence  time. 
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a.  ER  =  0.1 

Fig.  1  Enthalpy/Temperature  Map  for  H2/02 
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a.  Temperature  versus  Time  for  H2/02f  ER  =  0.1,  Enthalpy  Constant 
Fig.  2  Comparison  of  Plug  Flow  Reactor  and  Perfect  Stirred  Reactor 
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TABLE  I 

HYDROGEN/AIR  REACTION  SET 


H2  +  O  =  H  +  OH 
H2  +  OH  =  H  +  H20 
02  +  H  =  O  +  OH 
H20  +  O  =  OH  +  OH 
H  +  H  +  M  =  H2+M 
O  +  O  +  M  =  02  +  M 
H  +  O  +  M  =  OH  +  M 
H  +  OH  +  M  =  H20  +  M 
Reaction  scheme  and  rate 


f  =  2.25  x  1012  e-3900/T 
f  =  2.20  x  1014  e-5138/T 
f  =  1.90  x  1014  e-90O7/T 
f  =  8.30  x  1013  e_9108/T 
f  =  5.0  x  1018  x-1.15 
f  =  4.7  x  1015  T--28 
f  =  5.3  x  1015  e140°/T 
f  =  1.2  x  10i?  e252/T 
data  taken  from  Ref.  12 
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TABLE  II 

METHANE/AIR  REACTION  SET 


<  1 

CH4  +  0  = 

CH3  +  OH 

f  = 

1.7  x  1013  e"4380/1. 

2 

CH4  +  H  = 

«  ch3  +  h2 

f  = 

6.3  x  1013  e-635°/T 

3 

CH4  +  OH 

=  ch3  +  h2o 

f  = 

2.8  X  1013  e-2500/T 

4 

CH3  +  02 

=  CHO  +  H20 

f  = 

101 1 

5 

CHO  +  OH 

=  CO  +  H20 

f  = 

5.0  x  IOI3 

6 

00  +  OH  = 

=  C02  +  H 

f  = 

3.1  x  101 1  e-3oo/T 

7 

O  +  H2  = 

OH  +  H 

f  = 

3.0  X  1013  e-4T60IT 

8 

H  +  02  = 

OH  +  0 

f  = 

2.24  x  1034  e844/T 

9 

O  +  H20  = 

=  OH  +  OH 

f  = 

8.44  x  1013  e9AT0l T 

10 

H2  +  OH  = 

=  H  +  H20 

f  = 

2.19  X  1013  e-2S90/T 

11 

H  +  OH  + 

M  =  H20  +  M 

f  = 

2.0  x  101 9  T-i 

12 

CHO  +  M  = 

=  H  +  CO  +  M 

b  = 

1.0  x  1013  e-i4>40°/T 

13 

ch4  +  m  = 

:  ch3  +  h  +  m 

b  = 

1.5  x  lO39  e-5°.300/T 

14 

n2  +  o2  = 

NO  +  NO 

f  = 

1.36  x  IOI4  e-37.750/T 

Reaction  scheme  and  rate  data  from  Ref.  19  except  for  the  last 
reaction  whose  rate  constant  is  taken  from  Ref.  16  for  the 
reaction  O  +  N2  =  NO  +  NO.  * 
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TABLE  III 

SENSITIVITY  COEFFICIENTS  FOR  THE  METHANE/AIR  SYSTEM 


t.(d T/df,)  aK 


2100  2.5  2.3  4.2  -5. 
1900  1.0  0.8  2.8  -3. 


1700 

1500  0.7  0  5.8  9. 

2100  1.2  1.1  6.4  15.i 

1900  0.7  0.6  7.6  21.: 

1700  2.0  0.7  12.0  34.: 

1500  10.6  8.7  36.2  103. 

2100  1.8  1.3  16.0  81. 

1900  2.4  1.7  22.4  94.  i 

1700  7.1  5.4  45.3  149. 

1500  112.7  111.1  312.0  618.: 


0.8  2.8  -3. 

I  I 

.5  0*  3.2 


\  ~  <  |o.5 


2100  1. 
1900  2. 


5.2  21.3 

-2.0  36.0 

8.8  34.8 

5.1  44.7 

65.5 


58.9 

73.5 


7 

8 

2.2 

2.3 

0.8 

0.8 

0 

0 

0 

1.4 

0.8 

1.0 

0 

1.2 

0.7 

3.6 

12.6 

36.7 

2.9 

7.8 

10.5 

30.3 

202.9 

351.1 

0.8 


-3.6 
17.4 
6.6 
—9.7 
-68.0 
1 .3 
8.8  -18.2 
52.0  -57.8 


266.8  202.9  351.1  250.7  221.7  208.0  471.3-232.7 


-0.5 

0 

0 

-0.6 

0 

0 

-5.8 

0 

O 

-1.2 

-12.2 
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APPENDIX  III 

COMPUTER  PROGRAM  LISTING  AND  SAMPLE  OUTPUT 
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C  MAIN  PROGRAM  STIRRED  REACTOR  THEORY  I « T«  OSGERB Y 

REAL  MOOT V 
DIMENSION  A  (240) 

DIMENSION  X0(20 ) , X( 20 ) >  W(20 ) ,B(20 ) , BO (20 ) ,HM(20 ) , Y (20 ) , DCP(20) »BD 
•(20),CPM(20),Y0(20) 

COMMON  /SEC0ND/NN.N.N1.NR.NT, 1 T » M»L» N2 .KOUNT . NT1 

COMMON/RGAS/RG.RR 

COMMON  /EQL/F(20)  ,DM20>- ' 

COMMON/) IFTH/NNT 
COMMON/HRCON/ICR.NR2.NR3.NR4 
C  I  T  =  1  TEMPERATURE  CONSTANT 
ICR=80 
K=  0 
JJ=2 

100  CONTINUE 

REAU(29.3)  P.T.TO.DELT, THIN 
IK (P.EQ.O.O)  GO  TO  iOl 

READ( 29.2)  MDOTV.DMUV.NI .N.LT.MT, IT. NR,  NN.NNT 
READ (29, 1) ( XO ( I ), I  si. NT) 

READ ( 29, 1 ) ( X(  I), 1=1, NT  ) 

REAU(29,i)(h(I),I-l,NT) 

C  COMPUTE  SPECIES 
N1=N*1 
N2=N*2 
NT1=NT  *1 
NX =N 
MY=N*1 

20  RR=1. 98/26 
RU  =  82 . 0^67 
T  T  =  T 

XX  =  0.0 
YY=0.0 

C  NORMALISE  INPUTS 
DO  9  1=1. NT 
XX=XX+X(I) 

9  YY=YY+XO(  I  ) 

DO  10  1=1. NT 
X(I)=XU)/XX 
10  X0(  I)  =  XO(I )/YY 
HU  =  0  >  0 
Y Y=0 • 0 
W0  =  U • 0 
WM=0 . 0 
DO  4  1=1. NT 

C  CALCULATE  MOLEC.  NT. 

rfU  =  Wl)*XU(  I  )*W(  I  ) 

WW=HW+X( I )*H(  I  ) 

4  CONTINUE 

CALL  B93  ( TO , NT , HM , CPM, UCP , K ) 

C  CALCULATE  FtED  VALUES  OF  SPECIES  AND  ENTHALPY 

DO  iu  1=1, NT 
B0(1)=X0(I)/W0 
HO=HO+HM( I )*B0< I ) 

30  CONTINUE 

C  CALCULATE  FEED  DENSITY 

RO=P*WO/(8?.U56/*TO) 

WRITE  ( 21 , 3t> )  TO.  HU 
WRITE  (21,102) 

IF  (IT-1)  2tt.29.28 
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28  WRITE  <21,31) 

GO  T  3  14 

29  WHITE  <21,32) 

C  JfcFl-Nfc  MASS  FLOW  RATE  /  UNIT  VOLUME 
14  JO  12  L=1,LT 

II-  <<-.cQ.l)  GO  TO  13 
T=TT 

MU3TV=MU0TV*UMDV 
13  00  11  1=1, NT 

c  calculate  initial  reactor  species 

11  9<I)=X<1)/WK 

C  SET  TEMPERATURE  ESTIMATE 

DO  13  M=1.MT 
I) <M.EQ.1>  GO  TO  18 
T  =  T- CELT 

II-  <  (T-TMIN)  .LT.DELT  )  M  =  MT 
18  CAlL  833  <T,NT,HM,CPM,UCP.K> 

C  CHECK  TOR  NON-CONVEHGED'SOLUT ION  OR  BLOW  OUT 

II-  <  KO JNT .  EQ  •  60  )  GO  TO  22 
DO  21  J=N1,N2 
YU<J>=3< J)/HO< J) 

II-  <AbS<  Y0<J)-1.0)-d.UE-3)22.22.23 

21  CONT INJE 

22  UO  16  1=1, NT 
16  8  <  I)=X< I )/MH 

C  ESTIMATE  REACTOR  DENSITT  AND  RESIDENCE  TIKE 

23  R=3  /<B2.0567«T)*WW 

T  A J=  R/MUO  T  V 

CA.L  1  IMER  <  TT1 ) 

C A_L  CMMP  <H.80,P.T,R,Y.H,HC.CP,TAU#UU,XX.aU,A.NX,MY) 
CA_L  TIMiR  <  TT2 ) 

T IME=  TT2-TT1 

*R I T  5  <  21 «  34 )  TIME 

MR  I Te< 21.33)  N.NT. II , NR , KOUNT , M.L, NX , MY , MT , LT 

WHIIs  <21,39)  <U< I ) , 80  < 1 ),8D< I ),Y( I > , X  < 1) , XO  < 1 ) , 1=1, NT) 

WH I  I c  <  21 , 7 )  P,T,R.R-J.WM,HO,H,HO,CP 

WK  I  I E ( 21 ■ 17 )  TAU.MUOTV 

JT  =  1  T 


lh  <  NN )  43,15.43 

43  CA_L  RCHtK  <  T , B . bO ,H . H , T , H, HO , CP. T AU , N U , BU, A , NX, MY, MDOTV , XO . WO ,W ) 
IT  =  JT 


15  CONT INJE 
12  "ONT INJE 
GO  TM  107 
1  FORMAT (6cl2  ■  0  ) 
c  format < 2 Eio.o,bi5) 
3  FO4MAT<5-10.0) 

7  FORMAT ( y El2 • 4 ) 


17 

FORMAT  <2E12*4///) 

31 

FORMAT 

(1H0.3X.19H 

enthalpy 

CONSTANT) 

32 

format 

< 1H0 , 3X , 22H 

TEMPERATURE  CONSTANT) 

33 

FORMA  T (11 15 ) 

34 

FORMAT 

< 1E12.4) 

35 

FORMAT 

<3X,2E12.*> 

3V 

format 

(3X.6E12.4) 

102 

FORMAT 

< 1H1.3X.47H 

AU  I  At)  AT  I  C 

NELL-STIRRED 

101 

STOP 

END 

0000  ERRORS.  COMPILATION  COMPLETE. 


REACTOR 


CALCULATIONS) 


34 


AE  DC-TR-72-164 


SUBROUTINE  COMPia.BU.P,  I  ,R.  Y , H, MO , CP, T AU, WW , XX , BD , D , NX , MY > 
DIMENSION  B<s;0)»dQt20)»OCNX,HY),A<2Ui20>»CC20)>DBDB(2D»20),  DBDT  ( 20 
•).  D3DR ( 20 ) ,  BD(20 > . HMC20 > ,CPM< 20) , UCPTT20) , DP<20> , Y(20  )  .  Y0T20 

*>. VS(20> 

C0MM3N  /SECOND/NN.N.Nl.NR.NT.IT,  IUUMC3)  .KOUNT.HUUM(l) 

COMM3N/R3AS/RG.RR 

CA.L  d?l(AfC,BJ> 

<UN I = 1 
I  CT=  C? 

DE„T„=i) .  0 
DO  6  KC1.6U 
KO'JM  =  KN 

c  CALCULATE  partial  derivatives 

CA.L  B34(B. r,R,D6DB,UBDT,DBDR.BD) 

C  CA.C'JL ATE  MATRIX  COEF  . 

DO  1  J=1.N 
DO  1  1=1. N 
ZUsU.O 

DO  2  L=N1 , NT 

Zd  =  23+ ( DBD8 (J.L)-DHUR(J)*R*MH-BU(J)*HH)*A< I . L ) 

2  CONTINUE 

D( J. I)=-TAU*(DBUB(J, I )-UBDR< J)*R*WW-BD< J)*WH+ZB) 

1  CONTlNJt 
DO  4  1=1, N 
D(I,I)=1.0+U<1.1) 

4  D  ( 1 »MY)=30( I )-B  < 1 ) ♦ T AU*HD  <  I  ) 

C  SOi_VE  MATRIX 

CAlL  CHOLES  (D.NX.MY.U) 

26  IF ( NX )  12.13.13 

12  WH I TE ( 21 < 14 ) 

14  FORMAT (1H.3X.15HSI NUULAR  MATRIX) 

NX  =  N 

SO  T3  100 

13  DO  16  J*1 » N 
YO( J)=3( J) 

16  YS< J)=D( J.MY) 

C  TEST  FOR  CORRECTIONS  TOO  LARGE 

2PX=0.0 
DO  111  J  =  1,N 
ZPP=ABS(YS( J)/YU( J) ) 

IF  (2P3.GT.2PX>  2 PX=2PP 
7/7  CONTINUE 

2P  =  0 . 3  7/2PX 
IF ( 23 . GT . 1 . 0 )  2P=1.0 
C  CORHECT  SPECIES  CONCENTHATIONS 

DO  17  J=1,N 

15  YC J )  =  YO < J) ♦VSC J)*2P 
DO  17  1*1. N 
B(I)=Y(I> 

17  YCI)=(HCI)-BO(I))/( TAU*  8D(I))-1.0 
DO  10  J=  N1 , NT 

YU ( J)=3l J) 

B(J)SC(J) 

sot j)=o.a 

Y( J ) =0  .  U 
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JJ=J-1 
00  9  1=1. JJ 

BU( J)=dU( J)*A( I , J)*BU( 1) 

9  B( JJ =B( J)*A( I , J) *H( 1  ) 

IT  <dU<  J)  .tCJ.0.0)  U0  TO  667 
Y(  J)  =  n(  j)-bu(  J>  )/«  I  AU*BD  (  J  >  )  ”1  •  0 
66/  CONTlNJb 

IT  <dl]<J)>  o/,4b,4/ 

47  Y0<J)=d( 

GO  TO  3V 
4b  Y0(J)=1.J 

49  IT  (  d  (  J  1  )  4U.10.1U 

40  B(J)=1.U=-1U 
10  00 NT  I  rot: 

MM  =  (j .  Q 
UO  lb  1=1. NT 
lb  HW=K4  +  -j(I> 

4W  =  1.(J/N4 
R0_U=R 

r=3*ww/(RG*T ) 
taj=taj*r/rold 
C  ITsKATION  CONTROL 

UO  4 2  J=N1.NT 

IT  ( AH3< YO  C J)-l.Q)-1.0E-4>32.42.7 

41  CGNTINJb 
GO  TO  20 

1  IT  UP-1.0)  49.50.50 

50  IT  (11-11  32. 49. 52 

52  CAuL  dd4(T,NT , HH.CPN.UCHT . KUNT  ) 
HsQ.U 
CP=U.O 
DO  54  1=1, Ml 
ri=H*HM(  1  )*t3(  I ) 

54  0P=L,3*:PM(I)»B(1) 

C  TT-IPeM^lURfc  CONTROL 

IT  (ICT-0)  =4.55.54 

55  0b_f=(40-Hl/CP 

IT  ( UbLT*UfcLTL )  56.36,5/ 

57  T=T+UbLI*0.N4 
GO  f*  54 

56  T=r*UfcLT»U.4tt 
54  I CT=1 

ICT= ICT-l 
0b_ 1 _= jfcLT 

IT  (T-3'J;'0.'J)  3o.5b,>9 

59  T=5U0O.U 
GO  n  5U 

5b  IT  (T-J09.0)  61. 60. 60 

61  T=3U0.J 

C  TbHP  4ITHIN  1  DbG  K 

6 J  IT  ( AHS(HO-H)-CP)  49,4V, 6 
C  CH=L<  ■S.S.«.  EQUATION 

49  i)0  b  I  =1 ,  N 

IT  (  AH5  (  Y  (  I  )  )  “3  .  lib-4  )  8.6.6 

8  CONTlNJb 
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60  TO  20 
6  CONTINUE 

20  If ( KOUNT . EQ • 60  >  HRITE(21.99) 

XX=0.0 
H=0 . 0 
CP=0-0 

CALL  BB3(T»NT,HM,CPM.DCPT,KUNT) 

00  4  1=1. NT 
XX=XX+8( I )«MM 
I )*B(  I) 

3  CP=CP+CPM( I )*8(  I  ) 

100  RETURN 

99  FORMAT ( 27H  NO  CONVERGENCE  IN  60  STEPS) 
END 

0000  ERRORS.  COMPILATION  COMPLETE. 


SUBROUTINE  B81  (A.C.B) 

DIMENSION  AC20.20).B(20),C<20> 
DIMENSION  NS(20.5).SC(20.5),J(5) 

COMMON /SEC0ND/NN.N.N1. NR .NT . IDUMC6) 

DATA  KTP/O/ 

IF  ( KTP )  3.3.4 

3  NL=NT-N 

READ  (29.98)  ( J ( 1 ) . 1 =1 . NL) 

DO  1  1=1. NL 
JJ=J(I) 

READ  (29.99)  <  NS(K» I ) >SC(K. I )  ,K=1.JJ) 

1  WR  I  T  E  (  2i  ,  999  )  (  NS(K.  I)  ,§C(K.  I’)  ,K  =  1.JJ) 
999  FORMAT  (5( I10.E12.4) ) 

4  DO  6  1=1.20 

C( I )=0. 

DO  6  K=1 » 20 
6  AC  I . K )  =  0 . 

DO  5  1=1. NL 
I  I  =  N*  I 
C( I  I )  =  0 . 

JJ  =  J(  I  ) 

DO  2  K  =  1 >  J J 
L=NS( K .  I ) 

C( I  I )  =  C( 1 1 )*B(L)*SC(K.  I ) 

ACL. I  1 )=-SC(K,  I  ) 

2  A( I  I , 1 1 )  =0 . 0 

5  CONTINUE 

99  FORMAT  (  8(15,  F5.0)) 

98  FORMAT  (515) 

KT==1 

RETURN 

END 

0000  ERRORS.  COMPILATION  COMPLETE. 
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SU9R3UT INE  B33( T .NT . HH I . CPMI . UCPTMI . K) 

DIMENSION  HHI  ( 1  >  .CPU  l  (11  .DCPTMK1) 

DIMENSION  C06(l4),C(14.20) 

DIMENSION  CC(6).TT(5) 

COMMON  /tQL/M20),DM20) 

RG  =  1 • 98 /26 
TT.N  =  ALOG(T) 

TLN  =  l.-TTLN 
It  t  K )  4.7.4 
7  WRITE  (21.103) 

103  FORMAT  (1H1.3X.34H  I HERHO  DATA  CHECK  DATA  AT  TO) 

REAU  (29. 70)((C( I.KK). 1=1,14) ,KK=1, NT) 

/o  format (  /eh  .  o ) 

4  Tt  =  T 

IF  (Tb-300. I10.20.2u 
10  MR  I TE( 21 , 80 )  T 

80  FORMAT  (54H0  TEMPERATURE  IS  OUT  OF  RANGE  OF  CURVE  FIT  TABLES.  T 
HPel2  •  5) 

CA_L  EXIT 

20  I)  (Tb-lJOO. )30.40,40 
30  KL=0 

GO  TO  60 

40  It  (TE-5000.)  50,50.10 
50  KL=/ 

60  CONTINUE 
TT (1 ) =T 

00  1  M=2, 4 

1  TT (M) =  TT(M-1)*T 
TT(5)*1./T 
00  100  KK=1.NT 
00  5  J=l. 14 

5  C0£( J)*C( J.KK) 

CC(1)=  COECKL+l) 

DO  3  Js2, 6 
KKK=(KL«J) 

3  CC(J)=  COE(KKK)*TT( J-l) 

H=  CC(l)*.b*CU(2)*CC(3)/3.+.25*CC(4)*.2*CC(5)+CC(6) 

FX  = ( CC( 1 )*TLN  )-.5*CC(2)-(CC(3)*.5*CC(4) )/6.-.05*CC( 5>+CC(6)- 
*COs(KL*/) 

CP  =  (CC(1)+CC(2)+CC(3)*CC(4)+CC(5)) 

DC3=CC(2)+2.*CC(3)+3. *CC ( 4 ) +4 • *CC( 5 ) 

C  SPECIES  enthalpy 
HMI (KK)=H*RG*T 
C  SPECIES  HEAT  CAPACI f Y 

CPMl  ( KO=CP*RG 

C  SPECIES  HEAT  CAPACITY  DERIVATIVE 

DC3TMI(KK)=UCP*RG/T 
C  SPECIES  rREE  ENERGY 

F(<K)='X 

C  SPECIES  FREE  ENERGY  DERIVATIVE 

Dt ( KK ) s-H/T 

IMK.E3.0)  WRITE  (21. /I)  HHI  (KK>  .CPHI  (KK )  ,F ( KK ) . DF (KK > 

/I  FORMAT  (1H0.1P4E15.5/) 

100  CONTINUE 
K=1 

RETURN 

END 

0000  ERRORS.  COMPILATION  COMPLETE. 
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SU3H9UT 1 Nt  CHOLESt A.N.M.MATSYH) 
D I MfcNS I  ON  A  (  N  .  M } 

II-  (Atl.lJ.NE.O.O)  GO  TO  47 
DO  3/  J=2,N 

II-  (At  J.l)  .EO.O.O)  GO  TO  67 
H-wla  =  J 
GO  13  27. 

67  CONTINUE 

GO  T 3  54321 
27  DO  5/  K=1 » M 

TfcHPsAC IFL1P.K) 

A ( ITLI=,<)=A(1.K) 

At  1 . K  }  =  TEMP 
•57  CONTINJE 
47  DO  2  J=2.K 

Atl, J)=A(1, J)/A(l,l> 

2  CONI INJE 

-  GO  6  1=2. N 
DO  l  J  =  2» M 

IKtMATSYH.EQ.QlGO  TO  49 
II-  (1-J>69,66.67 
49  IT (J.GT.I)GO  TO  69 

68  K=J-1 
SUN=0.0 

JO  3  13=1. K 

SU4=SUM*4( I. I R 1  *  A  C  IH,  J) 

3  CONTINUE 

A< I . J ) = A ( 1 , JJ-SUH 
SO  T3  7 

69  4=1-1 
SU4=Q.O 

DO  4  13*1,  K 

SU4=SU4*A( I.IR)tA(IK.J) 

4  CONT I N jh 

IT (A((,|).EG.O.O)  GO  TO  54321 
A( I»J)  =  (A(  I , J ) -SUM  I  / A (  I  ,  I  ) 

GO  TO  J 

67  At  I , J ) = At  J, I )*A( J, J) 

7  CONTI N JE 
6  CONTINUE 
DO  52  <=2, N 
1=  4+i-< 

SU4=ti .  0 
LL  =  l  +1 

JO  51  I H=LL i N  1 
SU1=SU4*A(  I  ,  IR)*A(  IH.M) 

51  CONIlNJE 

A( 1,4)= A( I ,M)-SUH 

52  CONTINUE 

GO  13  12345 
54321  N=-l 
12345  RETURN 
END 

O'JOU  E9H3RS.  COMPILATION  COMPLETE . 
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5U3K3UT 1 NE  { 8 , T, K , DBDB, DBDT, DBDR , BD) 

DIMENSION  DBDB(20,2U) ,BU(20) ,UBUT(20).DBDR(20),B(20) 
COMM3N/SeCOND/NN,N,Nl,NH,NT,  I  T  ,  1 DUM( 3 ) .KOUNT , NT1 
C04M3N/R4TE/TT,FF{aU) ,BB(80) .OFF  (80)  .UBB(BO) 

DIMENSION  D(  30,18) ,  ll(8U,  7),ER  (BO,  7  >,KS(  80.7) ,  Ifc(80)  ,  IC(80) 
C0NM3N  /RRATE/A(80).EX(B0) , TEX (80) ,KK(8U,4), ITR(80>, ITPIBQ ), 
*  I B ( HU  ) 

DATA  KRU/07 
IF  (<HD)  1,1,2 
Nil  TtMP 
NT2  DENSITY 

NT  3  X 

1  NT1=NT*1 
NT  2  =  NT  *2 
N  T3*NT  *6 
KRi)  =  l 
B(N1 l)=l. 

RfcAU  (29,99)  NR 
00  10  IH=1,NR 
WRITE  (21,9/)  IR 

RfcAU  (29,99)  ITR(IR),1TP(IR),(KK(IR,I),1=1,4),A(IR),EX(1R), 
*TfcX(  IR) ,  It(IR)  ,  IC(  IK) ,  ItJ(  IR) 

WRITE  (21,999)  (KK(1R, I ) , I =1 , 4 ) , A ( I R ) , EX ( I R ) , TEX ( I R ) , 1 E ( I R ) . 
•ICC IR) ,IB( IK) 

999  FORMAT  ( 4 1 10 , 3E14 . 6, 4 1 5 > 

12  IF  ( IE( I R ) )  11.11.14 
14  L=IE( 1R> 

RfcAU  (29,96)  (Il(IR,I),ER(IR,l),I*l,L) 

WRIIc  (21.997)  (11(1R, I),fcR(IH,l), I=1,L) 

997  FORMAT  (3( I10.E12.4) ) 

11  00  16  1*1,4 

16  KS( (R. I )=KK( IR, 1 ) 

IF  (  1 C ( I R ) )  10.10,1/ 

17  M* I C  ( I R ) 

RfcAU  (29,97)  ( KS( IR. 1 ) , 1 =1 , M ) 

WRITE  (21,97  )  (KS( IR. I) , I=1.M) 

10  CONTI N JE 

2  CONTINUE 
TT  =  T 

as=o. 

DO  18  1=1, NT 

18  3S=BS*B( I ) 

IF  ( MOUNT . GT .1.AND.1T.EU.1)  GO  TO  866 
CA^L  HATES 
668  CONT IN Jt 

00  10U  I R=1 , NR 
JO  19  1*1, NT 

19  D(  IK  ,  l )  =  )  . 

K1=KK( I R , 1 ) 

K2  =  KM( IN. 2) 

K3*KMl IR.3) 
r<4  =  K<  (  IR.  4) 

IF  ( IT=( IR)-ITR( IR))  30,20,20 

20  IF  (ITR(IR)-l)  60,60.21 

C  EVALJATE  THIRD  BODY  EFFICIENCIES 
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60  CON  F  l  N  Jfc 
dSJM=H5 

If  ClbCIR))  0l.61|O2 

62  K=Ifc(  M> 

00  63  J=1,K 
JJJ=U(  I  rV  .  J) 

63  SS JH  =  BSUN  +  EH  C IK,J)»BCJJJ) 

61  l)l  =  f  f‘1  lK)«dCi<l> 

D2  =  -d«<  JR)*B<Ko) 

DCIH.Kl)=fF( IR)*bSUM 
DC IH.K3) =-8H< I R ) »dSUM 
DO  64  J= 1 , N I 

64  DC1R. J)=IUIH.J)+D1*U2 
If  (  I  b  (  1  R  )  )  c‘>.  65.66 

66  DO  6/  J=1,K 
JJJ=11( 13, J) 

67  J  C 1 R  »  J J J )  =D  ( lR,JJJ)+CUl+D2)*bR( 1R.J) 

65  DC IR.RT1 >=( Ul*D2)*BSUd 

DC  IR»NT2)=<  OFF  (  1R)*U1*DBBC  1R)*D2)*BSIIM 
DC1R,NT3)=J<  IR.NT1) 

*J0  T  3  10H 

C  EVALJAFb  TYPb  A*B=U*D 
21  CONriNJt 

DC  IH«K1)=FFC  [R)*b(K2> 

Dl-UC IR,Kl)*dCKl) 
DCIK.K2)=FF(IR)*BCK1)*D(IR,K2) 
U(1K,K3)=  -BBC IH)*B(K4) 

D2  =  U( 1R»K3)*BC<3) 

D  C I R • K4 )  =  -BBC IR)«3CK3)*DC IR.K4) 

J( IH.NT1)=D1*D2 

DC  I H » NF2 )  =  Dl*UFf C IR ) +U2*DB9 (IK) 

DC  I K ■ NT3 ) =DC  IR.NT1) 

30  13  10U 

C  EVALJATfc  TYPt  A*A*N=A2*M 

30  dSJM=bS 

If  ClbCIR))  311,311,31 

31  K= I t ( m 

JO  33  J=1 , K 
LL=UC  IK,  J) 

35  dSJM=B5UM*ERC IR, J)*BCLL) 

311  Jl  =  f  IR)*BCK1)*BCK2)*R 

J2=-dB( IR)»BCK3) 

JCIH.Kl)=f FC lR)*b(K2>*dSUN*R 
0CIrt.K2)=FF( IR)*B(Kl)*dSUM*H  ♦DCIR.K2) 
DC IR.K3) =-BBC 1R>  «BSUM 
DO  32  J=l,Nf 

32  DClH. J)  =  I>C  IK.  J)  +  D1*U2 
If  ClbCIR))  36,36,34 

34  DO  33  J=1,K 
LL= 1 1 C IR, J) 

33  DCIH,LL)=DCIR,LL)*CU1+D2)*ERCIR,J) 

36  DCIR.NT1)=CU1*U2)*BSUH 

DC IH,NT2)  =  CU1*DKF  C I W)*D2*UBB( 1 R ) )*9SUM 
DC IH.NT3>=C2.*D1*U2)*BSUM 
1U0  CONTlNJfc 
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C  CALCULATE  PARTIAL  DERIVATIVES 

DO  40  J=1.NT3 
DO  40  1  =  1, N 
40  DHDB 1 1 , J )  =  0 . 

DO  >0  IR=1 , NR 
K= I TR ( IR) 

L=  I  T3 ( I R ) ♦K 
M=K+1 

DO  50  J= 1 . NT3 
DO  51  1=1. K 
LL  =  KS( IR,  I  ) 

51  D8Dt)<LL,J)=UBDB(LL,J)-D(IR.J) 
DO  52  I=N,L 

LL=KS( IR. I ) 

52  DBDBILL. J)=DBDS(LL, J)*D( IR. J) 
i>0  CO VT  INUE 

DO  72  K=1,N 
DO  /I  JS1 , NT 

H  DBDHIK. J)=DBDB(K, J)»R 
BO(K)=DbOB(K,NTl)  *R 
DBDTCK)=DBDB(K,NT2)  *R 
DBDR(K)=DBDB(K.NT3) 

/2  CONTINUE 

V9  FORMAT  (  6 1  5  «  3E1 0 . 0  » «i  1 ) 

90  FORMAT  (3(  I&.E10.0)  ) 

97  FORMAT  (1015) 

RETURN 

END 

0000  ERRORS.  COMPILATION  COMPLETE. 
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SU3H0UT  I-NE  RATES 

COMMON /RATE/  T,FR(8U).BK(8Q>.UFRC80> .DBRC80) 

COMMON/RGAS/RR.R 
C0MM3N  /EQL/FC20) ,DF (20) 

COMMON  /RRATE/ AC  80 ) »  b < HO  )  .TEX  ( 80 ) , KK  ( 80 , 4 ) ,  ITR <  80  ) .  ITP (  80 )  , 
*18(80) 

COMMON/SEC0ND/NN,N,N1,NK,NT,  Il)UM(5).NTl 

C0MM3N/RRC0N/ICR,NR2.,NR.S.NR4 

RT=1 ./T 

RTT=  RT/T 

RRT=1 . / ( RR*T ) 

DRRT=-1./T 
F ( NT1 ) =Q • 

DF(NT1>=0. 

IF  (NN)  6./. 6 

6  AC ICR)=1.001*A(ICR) 

7  00  100  I R=1 , NR 
DdRC 1R)=0. 

UR ( I R ) =0  . 

DFR(IR)=0. 

F  R  (  I  R  )  =  A  (IR) 

IF  ( tr <  I R  )  )  1,2.1 

1  FR(IR)=FR(IR>*EXP(-b( IR)*RT> 

DFRClR)=c(IR)*RTT 

2  IF  (TEXCIR))  3,4.3 

3  FRl IR)=FR< IR)*T**TEX( 1R> 

DFRC IR>-0FR( IR>*TEX( IR)/T 

4  IF  (IBCIR))  5,100,5 

5  K=KK(1R,1) 

L=rfK(IR,2) 

J-KKC IR.3) 

I  =  KK( IR.4) 

EQ=F ( I>*F(J)-F(K)-FCLI 
IF  ( EQ* / 0 . )  100,9,9 
9  DE3=DF ( I ) *DF  C  J) -OF <K )-DF ( L  > 

EQ=bXPCEQ> 

IF  (IBCIR))  11,100,10 

10  0R(IR)=FG(IR>*EQ 
DBRC  IR)  =  DFR(  IRMDEQ 

IF  (ITPUR)-ITRC  IR))  8. 100. 100 

8  8RC IR)=dR( IR)*RRT 
DHRC IR ) =DBR( IR ) *DRRT 
GO  TO  100 

11  0=FR(IR)/<EQ*RRT) 

DH=UFR(IR)-0E0-DRRT 
BRC 1R)=FRC IR) 

DBR< IR ) =DFR( IR> 

DFR<IR)=DB 
FR{ I R) =0 
100  CONTINUE 

IF  (NN)  12.13.12 

12  AC ICR)=A( ICRJ/l.OQl 

13  RETURN 
ENO 

OOOO  ERRORS.  COMPILATION  COMPLETE. 


43 


AEDC-TR-72-164 


SUBROUTINE  RCHEK  ( T T . B, bO , P , R . Y , H, HO , CP. T A U, MH ,  BD.A.NX.MY.MDOTV 

•.xo.wo.vn 

REAL  MOOT V 
D  1  PENSION  A  (240) 

DIMENSION  XO ( 20 ) «  M(20 ) . B<20 ) .BU ( 20 ) ,HM(20 > , Y(20 ) ,DCP<20 ) »BD 

•(20).C=H(20).r0(20) ,DLE 1(80 >,8SAV(20). 20(20), DO (20) 

COMMON  /SEC0ND/NN»N*N1,NR»NT»IT»M»L*N2,K0UNT,NT1 

COMMON/)' I FTH/NNT 

COMMON/RGAS/RG.RR 

COMMON/RRCON/ ICR,NR2,NR3,NR4 

K=  1 

NR4=NR*1 
DO  20  1=1, NT1 
20 ( I )  =  X0 ( I ) 

D0( I >  =  30( I ) 

20  CONTINUE 

20  (NT1  )  =  Z0  (S>  ) 

DO  98  IHR  =  1»  NNT 
CALL  TIMER  (TTl) 

IF  ( I  MR-1 )  1.2,1 

C  INCREASE  SPECIES  NN 

1  Z0(NN)=2J(NN)*U.005«20(NT1) 

4TM=0.0 

DO  3  1=1, NT 

20(1)=2Q(I)/(1.0+20(NN)) 

HTM=MTM*20(I)*W(I) 

3  CONTINUE 

DO  4  1=1, NT 

DO ( I ) =20 (  I  )/WTM 

4  CONTINUE 

2  T  =  TT 

CALL  B33  (T,NT,HM,CPM,UCP,K) 

I  T*1 

R=®/(82. 0567«T)»WW 

TAU=R/MUOTV 

[CR=80 

C  CALCULATE  FEED  ENTHALPY 

CALL  COMP  (U,D0.P,T.R,Y.H,H0.CP,TAiJ,HM,XX,3D,A,NX,NY) 

IK  (K0JNT.E0.60)  GO  TO  39 
DO  19  1=1, NT 
GSAV ( I )  =  B(  I ) 

19  CONTINUE 
I  T  =  2 
HO=H 

C  CALCULATE  SENSITIVITY  COEF .  FOR  EACH  REACTION 

□0  t>  ICR*1, NR4 
IF  (KOUNT .E0.60)  GO  TO  6 
DO  7  J=N1,N2 

Y0( J)=S( J)/DO( J) 

IF  (ABS(YO(J)-1.0)-S».OE-3>  7,/,B 

7  CONTINUE 

6  DO  9  1=1, NT 

9  B(I)=aSAV(I) 

8  R  =  »  /<82.0567«T)*WH 

TAU=R/MUOTV 
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CALL  C3MP  (B.DO.P.T.R, Y. H , HO . CP . TAU, MW . XX . BD. A . NX , H Y ) 
OL=T( 1CP)=T 

It  { <OUNT . EO . 60 )  GO  TO  39 
DO  40  I =1 , NT 
3SAV( I )=S(  I  ) 

4U  CONT IN Jt 
5  CONT  1  N Jt 

3V  CONTlNJfc 

CA.L  TIMER  (TT2) 

TlMfc=TT2-TTl 
it  (IT-1)  23,29,23 

28  MR  I  T  e  (21,31) 

GO  T  3  33 

29"  MR  I  Tfc  (21,32) 

33  WRITE! 21, 11)  TIMfc 
DO  16  1  CR=1 , NR4 

16  DL=T( I JK)  =  (Di.ET<  ICR)-ULfcT(NR4) ) *1000.0 
C  SENSITIVITY  COtt.lS  DLE  t 

WR I Tfc ( 21 « 12 )  <  DLbT ( ICR) » ICR  =  1,NH4> 

MR  I  Tb (21, 13)  N , NT , I T , NR , KOUNT , M, L , NX , MY, I  MR 
WRITE  (21,14)-  (U(  1 )  ,D0(  1  ),BD(  1  ) ,  Y(  I ) ,  1=1, NT) 
wRITfc(21»12)P»T»R»RU»WW,WTM,H.H0,CP 
MR  I Tb ( 21 , 15  )  TAU.HDOTV 
ya  coNTiNjt 
RbTURN 

lj  FORMAT  (IH0.3X.19H  tNT HALPY  CONSTANT) 

11  FORMAT  (1E12-4/)  . 

12  FORMAT  ( 9512  •  4 ) 

13  FORMAT (915) 

14  FORMAT  (3X.4E14.6) 

15  FORMAT  (2E12-4///) 

31  FORMAT  (1H0.3X.1SH  ENTHALPY  CONSTANT) 

32  FORMAT  (1H0.3X.22H  ItMPfcRATURb  CONSTANT) 
tNJ 

0000  fcRRSRS.  COMPILATION  COMHLtTE. 
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DATA  CAHI3S 


■t* 

0\ 


Initial  Conditions 

p  r  ro  nixr  mis 

.2130.1)4  .3300.03  .2000.33  .5l)i)0*u3 

MIX1TV  DMI1V  NT  N  LT  MT  IT  NR  NN  NNT 

,1000.01  .2000.00  ts  •  1  1  1  14  0  0 

XCX1I  X(l)  Will 

.au; looio.oo  .20010000-04  .10000000.01 

.001100010*30  .11100000*03  .10000000*02 

.31'J03i)10*00  .2/763UOO-OJ  *1/300000.02 

,00300000*00  .32103000-02  .10310000*02 

.12*00030-02  .41100000-03  .10042000.02 

.00001000.00  .10000000-04  .13034000.01 

*30003100*00  .10200000-30  .20010000*02 

.0003,1300*00  .1O74OU00-03  .20010000*02 

.30001000*00  .17220000-03  .30000000*02 

.30030300*00  .13130000-02  .44010000*02 

.30301310*00  .43740000-04  .20100000*01 

.31323330-32  .34(00000-02  .32000000*02 

,1*«347U0-01  .20340000-01  .20010000*02 


TH=mo  D4T4  Coefficient.  ■ 

u  . 233o3UUU*01  *OOOOUUOU*UO  -30000000**0  *00000030*00 

**  .23301030*01  .00000031*00  ,00000000*00  .00010010*00 

,,  .307130*4,01  -.21737250-02  .37342203-05  -.2*947200-03 

”  .231733*7**11  -.104221*0-04  -.00017*20-00  .3*033020-11 

.34234/00-01  -.1110/230-02  .12466019-03  -.21333*30-0* 

u"  ,240*4344.01  .*<033061-03  -.21073*00-06  .1*002700-10 

„  _  .41363010*01  "-:i/2443O0-O2  .56902J20-03  -.43*30000-00 

.26  70  73  32  *  01  .3  g3l7l2q-02  -.0333l6o||-o0  .ll/*o0*o-0* 

ChT  .‘24*767**01  -.6*126600-02  .31602130-04  -.2*715430-0/ 

'  ,11/93/44*01  .10*503*0-01  -.40622100-05  .713/0230-0* 

CH  .33393032*01  ,426/0350-02  .20332600-06  -.11540330-03 

3  ,2032/662*01  .62304300-02  -.220*1300-03  ,300*3130-09 

2-HU  ,3464/1156.01  -.53/0060  0-03  .6*334300-03  -.66*30*30-00 

.32/00241.01  ,351*5000-02  -.137/2*00-03  .24/31430-0* 

rn  .3/3/1332*01  - .21/0*»30-02  ,30/57340-05  -.34/37/30-0S 

,2*51131**01  .155253/3-02  -.61911400-06  .11350300-09 

~ND  >4146*476*01  -.411*7230-02  .96*22467-05  -.70633640-00 

.31529360.1)1  .14339*55-02  -.5/3/04*0-0*  .10*70210-09 

rci  .21/01000*01  .10379l2o-01  -.10/33*00-04  .6345*100-03 

;  *  ,31129056*01  .31922*00-312.  -.A29/0200-05  .2414/403-0* 

H  .26461049*01  .41*3/120-02  -.9*119300-05  .V5122660-00 

2  .304310*7*01  .61167110-03  -.73*93300-00  -.20331*30-10 

.37109946*01  -.25167290-02  .03037353-05  -.02**0720-00 

*  .35*71129*01  .70145603-03  -.22336670-06  .424*0160-10 

. 36*1 H4K*ol  -.13337*00-02  .20503100-03  -.97600330-0* 

”2  . 26149761*01  ,15976320-02  -.62566200-36  .11315050-0* 


■nODDOOlO-Dl 

.234764*7.05 

-.460010  ’6*0*  3<XI<T<1000 

.1*0  0000)0*00 

. 25470497*05 

-.46001096*00  1000<  T<  5000 

.9077/550-12 

.2*13/1*0.05 

.2646,10/6.01 

-.55743600-19 

.29230007*09 

•49 4679-2*01_ 

-.-»25465o0-13 

.39092767*04 

.5625 50*9*00 

- .30452*00-15 

.36011792*04 

.55597040*01 

.14233690-11 

-.30260600*09 

- .611), 2 ,0.00 

-.61973600-14 

. 661343*1*0 1 

.99103900-11 

- .10166632*09 

-.917549*1*0.1 

-.4/^90000-13 

-.99996627*04 

•1250^944.02 

.41700400-12 

.19649793*09 

.2/03/4.8.01 

-.29276000-13 

.19707492*09 

.508411/0*01 

■205/7000-11 

-.27799993*04 

.46971067*01 

-.16727000-13 

-.27820147*04 

.7304423^*01 

. 7/2l*»00-l2 

-.l43633o**OS 

•26335499.01 

-.77003000-14 

-•14231627*09 

.65314450*01 

.22309910-11 

.25694290*01 

-.73720000-14 

.96927040*04 

.69446465*01 

-.16261000-11 

*.40352602*05 

.10664366*02 

-•1O74J000-1J 

-.48944043*09 

-. 72179769  *0j 

-.33093400-11 

-.96729400*03 

-.1411/650*01 

.2499J79I-14 

-.69491000*03 

-.1646.339*01 

.77062100-11 

-.10976700*04 

.3*00 )7  J  4*01 

-.J3460200-14 

-.11927910*04 

.37492699*01 

-.99772200-13 

-.10620300*04 

. 228749io*C1 

-.76097100-14 

-.69017449*03 

.639070  '9 »01 
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H 

THERMO  DATA 
hi 

5. 21069*04 

CH€CKc  DATA 

4 . 968^5*00 

AT  TO 

Fl/RT 

7.36022*01 

hj/RT2 

-2.91339-01 

O 

5.  93663*04 

5.234/4*00 

8.05398*01 

-3.33046-01 

OH 

9.32571*03 

7.14211*00 

-6.46677*00 

-5.21416-02 

h2o 

-5.77831*04 

8.02760*00 

-1.19645*02 

3.23075-01 

ch4 

-1.78732*04 

8.53504*00 

-5.24020*01 

9.99323-02 

ch3 

3.35080*04 

9.28114*00 

3.28309*01 

-1.87349-01 

CHO 

-3.16905*03 

8.27126*00 

-3.23501*01 

1.77187-02 

CO 

-2.64028*04 

6.96560*00 

-6.80687*01 

1.47623-01 

NO 

2.16132*04 

7.13246*00 

1.08962*01 

-1.20843-01 

co2. 

-9.40354*04 

8.89425*00 

-1.83457*02 

5.25768-01 

h2 

1.27436*01 

6.89381*00 

-1.57046*01 

-7.12515-05 

o2 

1.29930*01 

7.02368*00 

-2.466o6*01 

-7,2646i-o5 

n2 

1.28851*01 

6.96133*00 

-2.30323*01 

-7.20428-05 

TO  HO 

.300000*03  -.311001*02 
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BB1  DATA  CARDS 


3 

/sc 

.1000*01 

6 

'1000*01 

7 

•1000*01 

s 

.1000*01 

10 

.10.10*01 

11 

.1000*01 

4 

.1000*01 

5 

.2000*01 

6 

.1500*01 

1 

.50  j0*00 

3 

.sooo*oo 

7 

.5000*00 

12 

.iooo*oi 

10 

.1000*01 

7 

.5000*00 

2 

.5000*00 

3 

.50 j0*00 

4 

.0000*00 

6 

•5000*00 

9 

•5000*05 

13 

.1000*01 

9 

.5000*00 

BB4 

DATA  CARDS 

KK0R.4) 

AOR) 

EX(!R) 

TEX(IR) 

IE(IK) 

(C(IR) 

I  BOR) 

5 

2 

6 

3 

.170000*14 

•  430000*04 

.000000*00 

0 

0 

1 

5 

1 

4 

11 

.430030*14 

,435000*04 

.000000*00 

0 

0 

1 

5 

3 

4 

4 

.200000*14 

.oooooo*oo 

0 

0 

1 

6 

12 

7 

4 

.100000*12 

.300000*00 

>000000*00 

0 

0 

1 

; 

3 

0 

4 

.300000*14 

•900000*00 

.000000*00 

0 

0 

1 

3 

8 

10 

1 

.310000*12 

.300000*03 

.000000*00 

0 

0 

1 

n 

2 

3 

i 

.300000*14 

.474000*0* 

.000000*00 

0 

0 

1 

12 

1 

2 

3 

.224000*15 

.0441100*04 

.800000*00 

0 

0 

1 

4 

2 

3 

3 

.844000*14 

.947000*0* 

.000000*00 

0 

0 

1 

11 

3 

1 

4 

.219000*14 

•259000*04 

.oooooo*oo 

0 

0 

1 

3 

1 

4 

14 

.200000*20 

.000000*00 

'.100000*01 

0 

0 

1 

1 

0 

7 

14 

■144000*05 

.500000*00 

0 

n  -i 

1 

6 

5 

14 

•150000*20 

.503000*05 

•000000*00 

0 

0  -1 

12 

13 

9 

9 

.136030*15 

. 377500*05 

•000000*00 

0 

0 

1 
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SAMPLE  OUTPUT 

AUU9*T|r  VEIL -S'  tnkb.i  KEACTCR  CALCULATIONS 

TteKP-WAT JRg  CONSTANT 
.  i6t)3*0?  Time,  sec 

N  NT  IT.  NR  KK  M  L  NX  MY  MT  LT 

9  1J  1  14  >  1  |  9  II  4  1 

H  .«/3-)4B(I)  .nOOu+OuBOlI). 4297-01  -  .4214-0  3  41)  .7507-03  X(I). 0000*00  XOII) 

O  • 99J1-J4  .00011*00  ,1778*OuBDU>.io19-J3  .3133-02  >0000*00 


OH  .2384-jJ 

■ooou*au 

.4311*00 

.4200-J4 

.7770-02 

.0000*00 

HjO  .3324-12 

.0000*09 

.6011*01 

.3229-17 

.9227-01 

.0000*00 

CH,  ‘  .4U7-A 

•  1 764- 

-.3183*01 

-.1917-0/ 

.1151-03 

.4976-01 

CH3  .1420-04 

.0009*00 

.2930-91 

-.1U41-J4 

.5043-0 J 

.0000*00 

CHO  .1241-34 

.0009*00 

.2269-03 

.1258-32 

.2878-05 

.9000*00 

CO  -174U-JJ 

.0099*09 

.3148*08 

.6600-07 

.5247-02 

.0000*00 

NO  .1913-33 

.0009*00 

■1631*90 

. J0/2-U7 

.4822-02 

.0000*00 

COj  .15/0-02 

,oooo*ou 

.2839*01 

-.2154-07 

.4354-01 

.0000*00 

Ho  .4040-34 

.0900*00 

. 7343-01 

. 4  824- jB 

.1281-02 

.0000*01 

O2  -SS20-32 

.705  7-02 

-.6398*01 

.2860-05 

.9744-01 

.1991*00 

N2  .2457-31 

.2662-01 

-.9157-01 

. J203-J7 

•  74 J1 *00 

.7511*01 

>  34  00*0 1 P  .2100*04  T  . 5530 “0  3  p  .389>-02p„  .2003*02  W  ,2821*02  W.  .2079*03 H  -.3110*02  H0  .3314*00  CD 
.5530-0 3  tr  .1000*01  XXCI),  W,  NT  °  °  P 


ENTHAL^T  CONSTANT 
»161/«tl3 

ax  Mn|>2»17*0l  .2295*01  .4175*U1  -.5697*01  .1033*02 

J  “1. 2223*01  .9909*01  -.5595*01  >1487*02  -.1769*01 

9  13  2  14  1  1  1  9  10 

1  IMR 

-237939-04  . 000000*00  .429738-01  -.409549-03 

.963129-04  .900099*09  >17/766*00  .460856-04 

.238367- j3  >000000*00  .4310/1*00  .940294-04 

.332395-02  .000099*00  .601122*01  -.205791-05 

.436/1-*- 15  •  1  76436— C2  -.3l«267*01  .31072/-Q6 

.162029-34  .000009*00  .293015-31  .221412-04 

.125115-06  .000099*00  .225026*03  .196291-02 

.1/404 J-33  • OOOOOU«OU  .3X4750*00  -.264277-05 

>  1 31263- *i3  .000000*00  .183084*00  .293359-03 

.156962-02  .000000*00  .2BJB59*01  .256725-06 

.4J6019-04  .000300*00  > /3  4249-01  .229764-04 

•J51 978-02  ./05746-02  -.639766*01  .712246-05 

.265719-01  .266226-01  -.9l5422"Ul  .293359-03 

.3430*01  .211*0*04  .5530-03  -.0000*00  .2803*02 

.5530-03  .1000*01 


ENTHALPY  CONSTANT 
>1671*03 

AT  U  «f-»ll401  .22*9*01  .4172*01  -.569J*01  .1032*02  .1397*02  .2232*01  .2253*01  .2456*01 

J  1.^217*01  .9869*01  -.5595*01  .1894*02  -.1765*01  -0000*00 

9  13  2  14  1  l  1  9  10 

2 


.23/403-34 

.900000*09 

-430258-01 

--4l7041-OJ% 

.979376-  *4 

.000090*90 

.17/209*00 

.477565-04 

,231M2«-'J3 

.000003*00 

.431167*00 

.951352-04 

.  3  J  J  523-  1? 

.090000*00 

.603382*01 

-.206671-05 

.43351*-  .'5 

.17641J-02 

•318262*01 

.J19611-06 

.14351.-  J* 

.  8  42 1)44- 04  NN 

.136195-01 

.494010-04 

.124730-  .4 

.00000^*00 

.226938-03 

.199852-02 

.1745*4-  13 

.  00909*1*00 

.316437*00 

-.264813-05 

.  130B84-?3 

.090000*00 

.182403*09 

.292574-03 

■15//12-8? 

.000090*00 

.285233*01 

.255519-06 

.4}/5lS-j4 

.090000*00 

.737728-01 

■23461C-04 

.350307-32 

.705651-02 

.642296*01 

.711912-05 

.265586-31 

.266199-01 

.912317-01 

.292576-03 

.3430*01  . 21 u0*04  .5529-03  -.0000*00  «28C2*02  .2021*02  .2066*03  .2066*03  .3315*00 

.5529-03  .1000*91 


.1392*02  .2239*01  .2257*01  .2461*01 

.0000*00 


-.0000*00  .2079*03  .2079*03  .3314*00 
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